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We consider a quantum system coupled to a dissipative 
background with many degrees of freedom using the Monte 
Carlo Wave Function method. Instead of dealing with a den- 
sity matrix which can be very high-dimensional, the method 
consists of integrating a stochastic Schrodinger equation with 
a non-hermitian damping term in the evolution operator, and 
with random quantum jumps. The method is applied to the 
diffusion of hydrogen on the Ni(lll) surface below 100 K. 
We show that the recent experimental diffusion data for this 
system can be understood through an interband activation 
process, followed by quantum tunnelling. 
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The study of a quantum system coupled to a large 
background reservoir that leads to thermal fluctuations 
and dissipation in the dynamical evolution of the system, 
is of central importance in such fields as quantum optics 
l|l| , electronic conduction in nano-structures |], and dif- 
fusion of light adatoms on surfaces || . The standard 
formalism for this problem is through the master equa- 
tion for the density matrix psit) of the system . How- 
ever, this approach is not practical for condensed matter 
systems such as a hydrogen adatom moving on a metal 
surface. In this case, the density matrix would have di- 
mension N'^, where N is the product of the number of 
sites considered on the surface and the number of vibra- 
tional states included at each site. Typically, N would be 
at least of the order of 10^ rendering a direct numerical 
solution of the master equation unfeasible. 

Recently, an alternative approach known as the Monte 
Carlo Wave Function (MCWF) has been developed 
and applied to solve these type of problems in the field 
of quantum optics. In the MCWF approach, the evolu- 
tion of a quantum state is described by a stochastic 
wave equation, in which the original adiabatic Hamilto- 
nian Hs is only a part of the evolution operator: 
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Here the effect of each operator C^, acting on the quan- 
tum system represents a collision with the reservoir de- 
grees of freedom that takes the system from one quan- 



tum state to another. The new Hamiltonian iJ is non- 
Hermitian, built from Hg with an imaginary part added 
to account for dissipation: 

ih 
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The stochastic nature of quantum evolution is described 
by the quantities /o and j^,. They are random numbers 
such that the mean value of fp, is related to the scattering 
probabilities 
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with {f^)^5pfj,, and (/o)=l - 5p, where 6p = Y^p^ gives 
the probability for coherent propagation under H. With 
this choice of dynamics, it can be shown |^] that the 
quantity a{t) obtained by averaging a{t) = '^[t)){^{t)\ 
over all possible outcomes at time t of the MCWF evo- 
lution equation, coincides with the density matrix ps{t) 
obtained from the solution of the so-called Lindblad form 
of the master equation p|: 
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E^^^iC+C^Ps + PsC+C^ 



The equality between a and ps holds at all times t, pro- 
vided that it holds at t — 0. The particular form of the 
collision operators chosen in Eq. (|l|) is the most general 
one that preserves the normalization and positive defi- 
nitencss of the corresponding ps{t)- 

It is the purpose of the present Letter to demonstrate 
how the MCWF method can be used to tackle important 
transport problems in condensed matter physics in cases 
where the number of degrees of freedom is large enough 
(N ^ 10*) to make the density matrix approach unfeasi- 
ble. We consider here the case of a light adatom moving 
on a metal surface under conditions where classical acti- 
vated hopping rate between potential wells is negligible 
compared with the corresponding tunnelling rate. At 
present, there does not exist a clear consensus on the de- 
tails of the crossover from the classical activated behavior 
to the quantum tunnelling regime. In the Field Emission 
Microscopy (FEM) study for Ni and W substrates 
and in the latest STM study for H/Cu(001) §, a sharp 
crossover from classical diffusion to very weak tempera- 
ture dependence of diffusion was observed at a tempera- 
ture in the range of 60— 100 K. However, the Quasielastic 
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Helium Atom Scattering study for H/Pt(lll) yields 
no crossover down to T w lOOK. For the H/Ni(lll) sys- 
tem, recent optical studies showed a crossover behav- 
ior from the classical regime to a second activated regime 
with a lower activation energy below T « lOOK. This is 
in contradiction with the FEM data on the same system, 
which showed a crossover to a temperature independent 
diffusion at low temperatures [Q. Thus, while there is 
strong evidence that diffusion proceeds through quantum 
tunnelling at low temperatures, the detailed mechanisms 
for hydrogen diffusion on different substrates are not yet 
understood. Previous theoretical works do suggest that 
the details of the crossover is sensitive to the shape of 
the adsorption potential and not just determined by the 
barrier alone [|| ||. 

We will apply here the MCWF method to study the 
dynamics of H/Ni(lll). The low temperature activated 
behavior with a barrier of about 90 meV has been at- 
tributed to small polaron type activated tunnelling 
In our view, this is a highly implausible explanation. 
First, the polaron activation energy for H/Cu(001) ||] 
was determined to be ~ 3meV, then the relaxation en- 
ergy due to the adatom for H/Ni(001) has been calcu- 
lated to be 2.72meV [Q, and in our recent calculations 
for H/Pt(lll) [0 we also find a relaxation energy of just 
a few meV; the polaron activation energy is a fraction of 
the relaxation energy . We will show instead that the 
data can be explained in terms of tunnelling from the 
first excited vibrational states of the H adatom. 

We construct a semi-empirical potential U (r) based on 
available data as follows. The lowest energy adsorption 
sites are assumed to be the fee sites forming a 2D trian- 
gular lattice {1} with a lattice constant a = 2.581 A 
(see Fig. 1). Also, the neighboring hep sites at a distance 
of s = 1.49 A are taken to be equal in energy ||ic| ] 
(this is also supported by a recent ab initio calculation 
||14| ). Second, we fix the barrier between the fee and hep 
sites close to the value of 196 meV found in experiments 
. We use the vibrational excitation energy of 94 meV 
known from jl^, U{r) is constructed from local- 

ized Gaussians at both the fee and hep sites and adjust 
the Gaussian parameters, obtaining a fitting with a band 
gap between the centers of the Aq © Ai and the Eq © Ei 
bands of A = 96 meV and a separation between the low- 
est band and the top of the barrier between fee and hep 
sites of 207 meV. 

The adiabatic Hamiltonian Hs for our model is charac- 
terized by Bloch states {|k, to)} with corresponding en- 
ergy {ek,m}- Here to is the band index and k the 2D 
wave vector. The center positions and the bandwidths 
for the first few bands are listed in Table I. The first two 
branches form ID representations {Aq and Ai) of the 
symmetry group of the 2D triangular lattice, while the 
next four form 2D representations (Eq and Ei). 

We describe the H adatom as a linear superposition of 
energy eigenstates: 



|*(i)) -E^k,™(i)|k,m), (5) 

m , k 

with X^mkl^k.mP = 1- The frictional coupling to the 
substrate through electronic and phononic excitations is 
modelled by a general collision operator (|l|) , through 
which we model both intraband and intcrband transitions. 
It is represented as 

Cmi,m2.q F.^/^^^ q ^ |k + q, TOl) (k, TO2 | , (6) 
k 

where F is a (yet unspecified) transition rate, and fi in Eq. 
(1) now becomes a multiple index with two band indices, 
fi—{ml, to2, q). Thus the probabilities for scattering dpf^ 
are given by 

Sp, = {^it)\C+C^\^{t))St = J2 l^k,™J'r,.'5i. (7) 

k 

An important feature of the model is that for the low 
energy bands of interest, Aq® Ai and Eq® Ei, the com- 
posite bandwidths are much smaller than the energy gap 
A separating them (see Table I). This means that we 
need to consider only two kinds of transitions: interband 
transitions between the bands in the two groups, and in- 
traband transitions within each group. Since we do not 
have microscopic expression for the scattering rates Fjntra 
and Fintor we make one further simplification that is Fjntra 
= Pinter — P- Below, wc wiU show that the magnitude of 
D is controlled by the parameter 7 = HT/Ae, where A^; 
is the width of the upper composite band defined above. 

In our numerical calculations, the substrate is repre- 
sented by a 2D hexagonal box consisting of 180 x 180 
unit cells, with fully periodic boundary conditions. The 
size of the system is chosen such that the H adatom does 
not spread outside the boundary during the observation 
time t. To calculate the spatial a/3 elements of the tracer 
diffusion coefficient of H, we used the expression 

Dapit) = lim -^{{xa - {xa)o){S:p - {xp)o)), (8) 

where x is the position operator. The average (...) in Eq. 
(^) represents both the quantum mechanical average in 
a given state as well as the ensemble average over dif- 
ferent initial states. Statistical averages to compute D 
were performed with 1500 — 6000 initial states, for time 
intervals containing up to 10^ collisions. With a code 
parallelized on 4 processors, one point on the Arrhenius 
plot takes 2 — 4hrs, depending on the collision rates. 

The symmetry of the lattice implies that the diffusion 
tensor Daf3 is diagonal. Fig. 2 shows the temperature 
dependence of D for 7 = 1,5, and 10 on an Arrhenius 
plot. There is clear activated behavior D cx e^^"/*^^"^, 
with an activation energy Ea — 98.1 ± 0.5 meV. This 
is in excellent agreement with the experimental data of 
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Cao et al. shown in Fig. 2 as well, in the temper- 
ature regime below 100 K where Ea — 105 meV. Obvi- 
ously, with the inclusion of only the lowest bands in the 
present calculation, we cannot account for the classical 
high temperature region above 100 K where Ea — 196 
meV (l^. We can give a good qualitative description of 
the quantum regime, though, where the numerical results 
above indicate that the observed Arrhenius behavior for 
D corresponds to activated quantum tunnelling. 

The result for the temperature dependence can be un- 
derstood from the values of the bandwidths listed in Ta- 
ble 1. The bandwidths of the {£'o,i?i} states are more 
than one order of magnitude larger than for the lower 
bands (the delocalization was observed also in a recent 
experiment |p^ ). Thus, diffusion proceeds mainly via 
a collisions excitation to the upper band, followed by 
tunnelling to neighboring sites and de-excitation to the 
lower bands again. It is the Bose-Einstein factor niuj) 
(fiLD — A), needed to ensure detailed balance in thermal 
equilibrium ||l^ , that leads naturally to the activated Ar- 
rhenius behavior with an activation energy close to the 
energy gap A. Although the Arrhenius behavior of D 
does not depend on the ratio 7, its absolute magnitude 
is best fit to the experimental data by choosing 7 w 10. 
This should be taken only as an effective ratio between 
tunnelling and scattering, because e.g. polaron effects 
p8| , p^ which lead to a broadening of the levels and a re- 
duction in the tunnelling rate have been left out in the 
present calculation. 

The MCWF methods gives insight into the quantum 
dynamics by allowing to follow the dynamics of wave 
packets in real space and time. In Fig. 1 we show two 
typical trajectories, tracing the evolution of (f ) for a wave 
packet. The larger length scale for the trajectory at llOK 
refiects the larger value of the diffusion coefficient, which 
is due to a higher excitation rate into the upper bands. 
The trajectory at 70 K has points where the particle is 
in the ground state for a longer time and, by comparison 
to the trajectory at 110 K, it has less coherent propa- 
gation intervals in the upper band. The other point to 
note is that there are coherent propagation regions with 
tunnelling through several sites before a de-excitation. 
This can be quantified by studying the tunnelling length 
distribution Pg. We define the tunnelling length £ as the 
distance travelled by a wave packet in the upper band 
before it suffers a collision. It is found that asymptot- 
ically Pi decreases exponentially with ^, while it obeys 
a Poisson-like distribution at small values oi t {£ < s). 
This is similar to the jump distribution in the classical 
regime Regarding the dependence of D on 7, we 

have done simulations at T = 70 K and T = 110 K in 
the range 0.1 < 7 < 10 and found that _D cx 7^^ in this 
range. This inverse power law dependence on 7 is sim- 
ilar to the dependence of D on the microscopic friction 
77 in the classical regime |2^,^. However, the influence 
of the geometrical factor on the dependence of the jump 



distribution on 7 seems rather different from the classi- 
cal case. The crossover of the dependence on 7 or 77 from 
the quantum to classical behavior is a subject worthy of 
further investigations. 

To summarize, we have demonstrated through a model 
study of H diffusion on Ni(lll) that the MCWF method 
is a powerful tool in the study of quantum transport 
problems with many degrees of freedom. In addition, 
the real space nature of the method allows one to ex- 
tract interesting information about the dynamics of wave 
functions, not easily available in other means. As op- 
posed to the small polaron mechanism suggested earlier 
|l0| , our results suggest that the low temperature dif- 
fusion behavior observed in the work of Cao et al. for 
H/Ni(lll) jl^] has its origin in the tunnelling of the hy- 
drogen adatom from the first vibrational excited state. 
We plan to apply the same MCWF formalism to investi- 
gate other quantum diffusion systems, such as H/Pt(lll) 
§ and H/Cu(001) §], which show qualitatively different 
behaviors from H/Ni(lll) ||l^. The key is to start with 
a reliable adsorption potential through a combination of 
first-principle calculation and empirical inputs. 
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Bandwidths Ae™ and 


band centers for 


branches 1 — 


6. Groups 1 — 2 and 3 


— 6 form the compos- 


ite bands Aq < 
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m 
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2 (Al) 


0.008 


104.497 


3 {Eo) 


0.017 


200.346 


6 


0.017 


200.721 


4 (£1) 


0.146 


200.446 


5 


0.146 


200.621 



FIG. 1. (a) Trajectories at T = 70 K (smaller set) and 110 
K (larger set). 7 = 1 and the observation time was 3.1 x 10~^s. 
(b) Details of the path at T = 70 K. The black circles are 
excitations or de-excitations. Between two such consecutive 
points there are usually several random changes of the mo- 
mentum 

FIG. 2. Temperature dependence of D between 80 K and 
140 K, for 7 = 1,5, 10. The Arrhenius behavior is evident. 
The experimental data of Cao et al. [9] are shown for compar- 
ison. For 7 = 10, a prefactor Do of 2.71 x 10^ A /s is obtained. 
The experimental value of the prefactor Do is 2.4 x 10^ A / s 
[9] 
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